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A novel lattice approach is presented for studying systems comprising a large number of interact- 
ing nonrelativistic fermions. The construction is ideally suited for numerical study of fermions 
near unitarity-a strongly coupled regime corresponding to the two-particle s-wave scattering 
phase shift 8q = k/2. Such systems may be achieved experimentally with trapped atoms, and 
provide a starting point for an effective field theory description of nuclear physics. We discuss 
the construction of our lattice theory, which allows us to study systems of up to (but by no means 
limited to) 38 fermions with high accuracy and modest computational resources, and offer an 
overview of several applications of the technique. A more detailed discussion of applications and 
simulation results will be described in companion proceedings by A. N. N. and J-W. L. 
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1. Introduction 

Simulating theories at finite fermion density has been a long-standing challenge in lattice field 
theory. One of the main issues faced in Monte Carlo studies of such theories is the fermion sign 
problem: when a chemical potential is introduced, the action often becomes complex. A conse- 
quence of this is that the exponential of the action cannot be interpreted as a probability measure as 
required by standard Monte Carlo algorithms. In the case of lattice QCD at finite baryon number 
density, the sign problem has severely limited the explorable regions of the QCD phase diagram in 
the temperature-chemical potential (T-/j.b) plane to the regime where {J-b/T < 1. 

Alternatively, at zero temperature, finite densities can be achieved in a canonical ensemble 
setting by considering multi-fermion correlation functions in a finite box. In this case, however, a 
closely related problem emerges: correlators involving fermions typically have a signal/noise which 
decays exponentially with the time separation of the source and sink. This makes identification 
of effective mass plateaus difficult— if not impossible-for large numbers of fermions. In QCD, 
standard Lepage arguments ^ suggest that the signal/noise for multi-baryon correlation functions 
decay with a time constant T^ 1 ~ B(m p — 2>/2m n ), where B is the baryon number, m p the proton 
mass and m n the pion mass. 

In an effort to better understand the sign and signal/noise problems as well as their interrela- 
tionship, we choose to study multi-fermion systems in a much simpler setting than QCD. One of 
the simplest nontrivial and interesting theories describes a dilute two-component system of nonrel- 
ativistic fermions in the unitary regime. This regime corresponds to a two particle s-wave scattering 
phase shift near 8q = n/2, and when this limit is achieved, the bound G/ = o < 4n/p 2 on the s-wave 
scattering cross-section becomes saturated. Although these systems are dilute, they are strongly 
interacting and therefore require a nonperturbative treatment for reliable study. Unitary fermions 
are interesting not only because they can be studied experimentally using trapped ultra-cold atoms, 
but also because they serve as a starting point for a lattice effective field theory (EFT) descrip- 
tion of nuclear physics. The latter is due to the fact that the and Si scattering lengths for 
nucleon-nucleon scattering are unnaturally large compared to the range of interaction. 

In this work, we focus on a new canonical ensemble approach for simulating large numbers of 
nonrelativistic fermions in the unitary regime and at zero temperature. Using this new approach, 
we consider two systems in particular: unitary fermions in a finite box and unitary fermions in 
a harmonic trap. From numerical simulations of these systems, one may extract experimentally 
measurable non-perturbative quantities such as the Bertsch parameter (i.e., the ratio of the energy 
to that of the free gas energy in the thermodynamic limit) and pairing gap. Furthermore, moving 
away from unitarity one may test a set of universal relations involving a quantity known as the 
"integrated contact density", first discovered by Tan [Q, [| |J]. 

These proceedings focus primarily on the details of our lattice construction, as well as the 
simulation and parameter tuning methods used in our studies of unitary fermions. Many past 
numerical simulations at zero temperature have been variational in nature, offering only an upper 
bound on the Bertsch parameter and possessing unknown systematic errors on other quantities 
[||, ^, Our approach is nonvariational and is therefore in principle free from such systematic 
errors. However, as is common among all simulations of this type, obtaining reliable results for 
the spectrum and matrix elements requires a good choice of interpolating operators which possess 



2 



A new approach for studying large numbers offermions in the unitary regime 



Michael G. Endres 



large overlap with the states of interest. Details of how we construct optimal sources and sinks as 
well as consideration of systematic errors due to finite volume and lattice spacing effects are the 
focus of our companion proceedings [§, ^]. In those proceedings we also present results of our 
initial exploratory studies of up to 20 trapped and 38 untrapped fermions, and present preliminary 
values for the Bertch parameter and pairing gap. 



2. Lattice construction 



The starting point for our construction is a highly improved variant of the nonrelativistic lattice 
action first proposed in []To|], given by: 



(2.1) 



This action describes two species of interacting fermions y = defined on a T x L 3 lattice 

with open boundary conditions in the time direction and periodic boundary conditions in the space 
directions. The derivative d t represents a forward difference operator in time and V 2 is a non-local 
lattice gradient operator which we define so-as to give a "perfect" continuum-like single particle 
dispersion relation for free fermions. A four-fermion contact interaction is achieved via a Gaussian 
or Z2 auxiliary field associated with the time-like links of the lattice. The operator C acts only in 
space, and in principle may include derivative interactions to an arbitrary even order in momenta. 
We may express Eq. 2.1 succinctly as S = yKy, where the time components of the fermion 



matrix K are given in block-matrix form by: 
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with 



(2.3) 



Note that the L 3 x L 3 matrices D, X and C act only in space and that 0(t) is a diagonal matrix with 
independent random elements. In momentum space the specific expressions we use for D and C 
are: 



<P|0|P'> 



*P 2 /(2M) 



P ' IpI < a 
IpI > A 



(p|c|p / ) = c( P )V 



(2.4) 



where A = 71 is a hard momentum cutoff imposed on the fermions and C(p) is some analytic 
function of p 2 which may be determined order by order in momenta from scattering data (see 
Sec. H for details). 
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Fermion propagators from time slice zero to time slices x may be expressed as a sequence of 
applications of D and X operators: 

K-\x;0) =D- l X{x-l)D- 1 ...D- l X(0)D- 1 , (2.5) 

which provides a simple recursive approach for computation. Inversion of the non-local D operator 
and application of X(x) may be performed efficiently with fast Fourier transforms (FFTs); it is this 
feature that allows us to use a perfect dispersion relation and momentum dependent interactions. 



In the free-field limit, one may explicitly verify that our definition for D provided in Eq. yields 



a perfect dispersion relation for fermions. In momentum space the fermion propagator reduces to: 

(p\K-\x;0)\p') =^ £(P)T V , E(p) = i , (2.6) 
for momenta less than the cut-off A. 



3. Simulation method and cost 

Since K is an upper tri-diagonal block matrix, the fermion determinant obtained by "integrat- 
ing out" the fermions is given by detK = detD r . We see that the determinant is independent of the 



auxiliary field, and therefore the full numerical simulation of Eq. 2. 1 is equivalent to a quenched 



simulation. By construction, our simulations are therefore free of the sign problem. Furthermore, 



the decomposition of fermion propagators as a product of random matrices in Eq. ^5] provides an 
unusual interpretation for our approach: in essence we perform Euclidean time-evolution of sin- 
gle particle wave functions over random background noise (either or Gaussian). Multi-particle 
sources and sinks may be constructed from a direct product of single particle wavef unctions, and 
in the case of sinks, more elaborate constructions may be considered as well which incorporate 
pairing correlations [||, . 



Numerical simulation of Eq. |2JJ consists of four steps: 1) lattice generation, 2) propagator gen- 
eration, 3) projection of time-evolved states onto sinks and 4) computation of Slater determinants 
(i.e., anti-symmetrization of initial and final states). The relative computational cost of each of 
these steps on an L = 32 lattice is shown in Fig. |j] as a function of the number of identical fermions. 
To give meaning to the vertical axis, note that the 0(N°) curve indicates the numerical cost of 
generating ^(L 3 ) random numbers (i.e., lattice generation). The steps l)-3) scale like T x L 3 or 
T x L 3 logL 3 , whereas step 4) is independent of the spatial volume. Although step 4) scales like 
0(N 4 ), 1 the computational cost of this step is negligible even for as many as N ~ 50— 100 identical 
fermions on the typical volumes we consider, which range from L ~ 12 — 64. 

4. Transfer matrices 



We may translate our lattice action Eq. gT| into Hamiltonian language with a suitable reinter- 



pretation of the expression for multi-fermion correlation functions. Such correlation functions are 



1 The cost of computing determinants scales like JV 3 , however, we compute determinants of all N sub-matrices as 
well, giving rise to an additional power of N. 
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Figure 1: Cost of numerical simulations as a func- 
tion of the number of fermions N. 




Figure 2: A plot of the three-dimensional £ function 
S(rj) given by Eq. 5.2 



obtained from an ensemble average of direct products of propagators. Since the propagator defined 
in Eq. ^ is itself a product of uncorrected random matrices (the auxiliary field is action-less), the 
multi-fermion correlation function will factorize into a matrix product of ensemble averages. If we 
define: 

^ = ^ 1 /2( 1 _^)^l/2 j (41 ) 

where 

3i = D- l l 2 ®...®D- l l 2 , {l-V) = (X(t)(8)...(8)Z(t)) (4.2) 

' * ' * * ' 

N N 

are V N dimensional matrices, then the N-fermion correlator may be written as: 

(K-\t,0)®...®K- 1 {t,0)) =^-1/2^^-1/2 > (4 _3) 
* v ' 

N 

and we may identify EF as the transfer matrix and M' = — log 2? as the Hamiltonian of the N- 
fermion system. 

In the case of two fermions (one up and one down), the transfer matrix evaluated in momentum 
space is given by: 2 

(P q l^lpq) e _ {p 2 +p ,2 W+ql 2 )/m , C(p) - —l f C 2n 2n (p) , (4.4) 

where we have expanded C(p) in the operator basis: 

2 n{v)=M n (\-e-v 2 l M y ^\> 2n , for |p| << 1 . (4.5) 

This transfer matrix can be diagonalized exactly, with the zero center of momentum energy eigen- 
values (e~E) given by solutions to the integral equation: 

Z,C 2n I 2n (p) = 1 , I 2n (p) = 1 £ (4.6) 

n v qeBZ e 1 

2 We have recently implemented a new Hermitian, Galilean invariant and analytic version of this interaction which 
corresponds to replacing: ^/C(p')\/C(q') —¥ C(p — p'); results using this improved interaction will be presented in 
detail in a future publication. 
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where p = y/ME and q are momenta within the first Brillouin zone (BZ). 



5. Lattice parameter tuning and the continuum limit 



In the continuum, Luscher's formula allows us to relate the two particle scattering phase shifts 



at infinite volume to the discrete energies of the same system at finite volume [ ]11| , |12| ]. Having 
solved the two body problem exactly on the lattice at finite volume, we may now relate the lattice 
couplings of our theory to continuum scattering data. Starting from the effective range expansion 



p cot do 



1 1 2 

a 2 



(5.1) 



one can relate the scattering length (a), effective range (ro) and higher order shape parameters to 
the s-wave scattering phase shift. Given an expression for /;>cotSo one may then determine the 



continuum energy eigenvalues in a finite box from the solutions to [13] 



pcotSo 
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TlL 



S(t?), 5(tj) 
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|n|<A 
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1. 



(5.2) 



,k) defined 



with X] = (pL/2n) 2 . Finally one may tune the lattice couplings C^n (for n 
in Eq. B3I by matching the lattice eigenvalues predicted by Eq. 4.6 to the lowest k continuum 



energies predicted by Luscher's formula. In this way, we may absorb all temporal and spatial 
lattice discretization errors into our definition of the couplings. 

The continuum limit for our lattice theory corresponds to the limit of infinite scattering length 
as measured in units of the lattice spacing. In order to maintain a finite physical scattering length, 
one must take this limit while keeping other physical quantities measured in units of the scattering 
length held fixed. In the case of unitary fermions, however, we need not concern ourselves with 
such complications and simply take pcotdo = 0. We therefore tune the couplings C2 n so that the 
lattice eigenvalues match the lowest k roots of S(rj) shown in Fig. ||| 

In Fig. H we plot the percent deviation between the roots of S(rj) and the lattice eigenvalues 
predicted by Eq. ^1] for up to seven tuned C- values for unitary fermions. For eigenvalues less than 
the k-th, the deviation is exactly zero by construction. For eigenvalues beyond the k-th, we find 
that the percent deviation remains quite small even for eigenvalues as high as 27, corresponding to 
three filled shells. Given a tuned set of C-values, we may also take the eigenvalues predicted by 
Eq. 4^6 and insert them back into Luscher's formula, giving a lattice prediction for /?cot5o. Fig. U 
shows a plot of the predicted pcot^o for up to seven tuned C-values. We find /?cot5o <C 1 for a 
wide range of momenta, extending well beyond that of the k-th eigenvalue we tuned to. 



6. Conclusion 



We have developed a new approach for simulating a large number of nonrelativistic fermions 
in the unitary regime. In these proceedings we have described some of the details of our lattice 
construction, an efficient numerical implementation of the theory and a method for tuning the 
lattice couplings to scattering data. Application of these ideas to unitary fermions in a finite box 
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Figure 3: Percent deviation in 77 between exact lat- 
tice eigenvalues (using M = 5 and L = 32) and con- 
tinuum Luscher eigenvalues. 
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Figure 4: Implied pcot^o obtained from exact lat- 
tice eigenvalues (using M — 5 and L — 32) and 
Luscher's formula. 



[]|] and in a harmonic trap 3 |§] are discussed in greater detail in our companion proceedings. There 
the issue of finding optimal sinks/sources as well as finite volume and cutoff effects are explored, 
and preliminary results for the Bertsch parameter and pairing gap reported. 
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